#include "catch.hpp"
#include "translational/translational_util/sbfill.h"

TEST_CASE( "sbfill" ){
  GIVEN( "invalid inputs" ){
    int ndmax = 5100, nbt = 2545;
    double delta = 1.791435E-3, be = 0.15;
    std::vector<double> sb (ndmax,0), ap (ndmax,0), betan {0.15, 0.18, 0.22};
    for ( size_t i = 0; i < ap.size(); ++i ){ ap[i] = (i+1)*0.0001; }
    sbfill(sb, nbt, delta, be, ap, betan, ndmax );

    std::vector<double> correct {0, 0, 0, 0, 0, 2.99592E-4, 2.94201E-4, 
    2.88907E-4, 2.83708E-4, 2.78602E-4, 2.73589E-4, 2.68666E-4, 2.63831E-4, 
    2.59083E-4, 2.54421E-4, 2.49842E-4, 2.45347E-4, 2.40931E-4, 2.36596E-4, 
    2.32338E-4, 2.28157E-4, 2.24051E-4, 2.20020E-4, 2.16060E-4, 2.12172E-4, 
    2.08354E-4, 2.04605E-4, 2.00923E-4, 1.93916E-4, 1.86053E-4, 1.78509E-4, 
    1.71271E-4, 1.64327E-4, 1.57664E-4, 1.51271E-4, 1.45138E-4, 1.39253E-4, 
    1.33607E-4, 1.28190E-4, 1.22992E-4, 1.18005E-4, 1.13221E-4, 1.08630E-4, 
    1.04225E-4, 9.99999E-5, 9.59453E-5, 9.20551E-5, 8.83227E-5, 8.47415E-5, 
    8.13056E-5, 7.80090E-5, 7.48460E-5, 7.18113E-5, 6.88996E-5, 6.61060E-5, 
    6.34257E-5, 6.08540E-5, 5.83866E-5, 5.60193E-5, 5.37479E-5, 5.15687E-5, 
    4.94778E-5, 4.74716E-5, 4.55468E-5, 4.37001E-5, 4.19282E-5, 4.02282E-5, 
    3.85971E-5, 3.70321E-5, 3.55306E-5, 3.40900E-5, 3.27078E-5, 3.13816E-5, 
    3.01092E-5, 2.88884E-5, 2.77171E-5, 2.65933E-5, 2.55150E-5, 2.44805E-5, 
    2.34879E-5, 2.25356E-5, 2.16218E-5, 2.07451E-5, 1.99040E-5, 1.90970E-5, 
    1.83227E-5, 1.75798E-5, 1.68670E-5, 1.61831E-5, 1.55269E-5, 1.48974E-5, 
    1.42933E-5, 1.37138E-5, 1.31577E-5, 1.26242E-5, 1.21124E-5, 1.16213E-5, 
    1.11501E-5, 1.06980E-5, 1.02642E-5, 9.84809E-6, 9.44879E-6, 9.06567E-6, 
    8.69810E-6, 8.34542E-6, 8.00705E-6, 7.68240E-6, 7.37090E-6, 7.07204E-6, 
    6.78530E-6, 6.51018E-6, 6.24622E-6, 5.99296E-6, 5.74997E-6, 5.51683E-6, 
    5.29315E-6, 5.07853E-6, 4.87261E-6, 4.67505E-6, 4.48549E-6, 4.30363E-6, 
    4.12913E-6, 3.96171E-6, 3.80108E-6, 3.64696E-6, 3.49909E-6, 3.35722E-6, 
    3.22109E-6, 3.15837E-6, 3.28595E-6, 3.41868E-6, 3.55677E-6, 3.70045E-6, 
    3.84992E-6, 4.00544E-6, 4.16723E-6, 4.33556E-6, 4.51069E-6, 4.69290E-6, 
    4.88246E-6, 5.07969E-6, 5.28488E-6, 5.49835E-6, 5.72045E-6, 5.95153E-6, 
    6.19193E-6, 6.44205E-6, 6.70227E-6, 6.97300E-6, 7.25467E-6, 7.54772E-6, 
    7.85260E-6, 8.16980E-6, 8.49981E-6, 8.84315E-6, 9.20036E-6, 9.57200E-6, 
    9.95865E-6, 1.03609E-5, 1.07794E-5, 1.12148E-5, 1.16678E-5, 1.21391E-5, 
    1.26295E-5, 1.31397E-5, 1.36704E-5, 1.42226E-5, 1.47971E-5, 1.53949E-5, 
    1.60167E-5, 1.66637E-5, 1.73368E-5, 1.80371E-5, 1.87657E-5, 1.95237E-5, 
    2.03124E-5, 2.11329E-5, 2.19865E-5, 2.28747E-5, 2.37987E-5, 2.47600E-5, 
    2.57601E-5, 2.68007E-5, 2.78833E-5, 2.90096E-5, 3.01814E-5, 3.14006E-5, 
    3.26690E-5, 3.39886E-5, 3.53615E-5, 3.67899E-5, 3.82760E-5, 3.98222E-5, 
    4.14307E-5, 4.31043E-5, 4.48455E-5, 4.66569E-5, 4.85416E-5, 5.05024E-5, 
    5.25424E-5, 5.46648E-5, 5.68729E-5, 5.91703E-5, 6.15604E-5, 6.40471E-5, 
    6.66342E-5, 6.93258E-5, 7.21261E-5, 7.50396E-5, 7.80708E-5, 8.12244E-5, 
    8.45053E-5, 8.79189E-5, 9.14703E-5, 9.51651E-5, 9.90092E-5, 1.03008E-4, 
    1.07169E-4, 1.11498E-4, 1.16002E-4, 1.20688E-4, 1.25563E-4, 1.30635E-4, 
    1.35912E-4, 1.41402E-4, 1.47114E-4, 1.53056E-4, 1.59239E-4, 1.65671E-4, 
    1.69228E-4, 1.72021E-4, 1.74859E-4, 1.77745E-4, 1.80678E-4, 1.83660E-4, 
    1.86691E-4, 1.89772E-4, 1.92903E-4, 1.96087E-4, 1.99323E-4, 2.02612E-4, 
    2.05955E-4, 2.09354E-4, 2.12809E-4, 2.16321E-4, 2.19891E-4, 2.23519E-4, 
    2.27208E-4, 2.30957E-4, 2.34769E-4, 2.38643E-4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

    for (int i = 2500; i < 2760; ++i){
      REQUIRE( correct[i-2500] == Approx(sb[i]).epsilon(1e-6) );
    }
    for (int i = 0; i < 2501; ++i){
      REQUIRE( 0.0 == Approx(sb[i]).epsilon(1e-6) );
    }

    for (size_t i = 0; i < ap.size(); ++i){ ap[i] = (i+1)*0.0001; }
    for (size_t i = 0; i < sb.size(); ++i){ sb[i] = 0.0;  }

    be = 0.015;
    betan = {0.015, 0.18, 0.22};
    correct = {0.0, 2.976482E-4, 2.922919E-4, 2.870320E-4, 2.818668E-4, 
    2.767946E-4, 2.718136E-4, 2.669223E-4, 2.621189E-4, 2.574020E-4, 
    2.527700E-4, 2.482214E-4, 2.437546E-4, 2.393682E-4, 2.350607E-4, 
    2.308307E-4, 2.266768E-4, 2.225977E-4, 2.185920E-4, 2.146584E-4, 
    2.107956E-4, 2.070023E-4, 2.032772E-4, 1.998421E-4, 1.983438E-4, 
    1.968567E-4, 1.953808E-4, 1.939160E-4, 1.924621E-4, 1.910191E-4, 
    1.895870E-4, 1.881656E-4, 1.867548E-4, 1.853547E-4, 1.839650E-4, 
    1.825857E-4, 1.812168E-4, 1.798582E-4, 1.785097E-4, 1.771714E-4, 
    1.758430E-4, 1.745247E-4, 1.732162E-4, 1.719175E-4, 1.706286E-4, 
    1.693493E-4, 1.680797E-4, 1.668195E-4, 1.655688E-4, 1.643275E-4, 
    1.630954E-4, 1.618726E-4, 1.606590E-4, 1.594545E-4, 1.582590E-4, 
    1.570725E-4, 1.558948E-4, 1.547260E-4, 1.535660E-4, 1.524147E-4, 
    1.512719E-4, 1.501378E-4, 1.490122E-4, 1.478950E-4, 1.467861E-4, 
    1.456856E-4, 1.445934E-4, 1.435093E-4, 1.424334E-4, 1.413655E-4, 
    1.403056E-4, 1.392537E-4, 1.382096E-4, 1.371734E-4, 1.361450E-4, 
    1.351243E-4, 1.341112E-4, 1.331057E-4, 1.321078E-4, 1.311173E-4, 
    1.301343E-4, 1.291586E-4, 1.281902E-4, 1.272291E-4, 1.262753E-4, 
    1.253285E-4, 1.243889E-4, 1.234563E-4, 1.225307E-4, 1.216120E-4, 
    1.207003E-4, 1.197953E-4, 1.188972E-4, 1.180058E-4, 1.171210E-4, 
    1.162429E-4, 1.153714E-4, 1.145064E-4, 1.136479E-4, 1.127959E-4, 
    1.119502E-4, 1.111109E-4, 1.102778E-4, 1.094510E-4, 1.086304E-4, 
    1.078160E-4, 1.070076E-4, 1.062054E-4, 1.054091E-4, 1.046188E-4, 
    1.038345E-4, 1.030560E-4, 1.022833E-4, 1.015165E-4, 1.007554E-4, 
    9.999999E-5, 9.925025E-5, 9.850614E-5, 9.776760E-5, 9.703460E-5, 
    9.630709E-5, 9.558504E-5, 9.486841E-5, 9.415714E-5, 9.423117E-5, 
    9.477307E-5, 9.531807E-5, 9.586622E-5, 9.641751E-5, 9.697197E-5, 
    9.752963E-5, 9.809049E-5, 9.865457E-5, 9.922190E-5, 9.979249E-5, 
    1.003663E-4, 1.009435E-4, 1.015240E-4, 1.021078E-4, 1.026950E-4, 
    1.032856E-4, 1.038795E-4, 1.044769E-4, 1.050777E-4, 1.056820E-4, 
    1.062897E-4, 1.069010E-4, 1.075157E-4, 1.081340E-4, 1.087558E-4, 
    1.093813E-4, 1.100103E-4, 1.106429E-4, 1.112792E-4, 1.119191E-4, 
    1.125627E-4, 1.132100E-4, 1.138611E-4, 1.145158E-4, 1.151744E-4, 
    1.158367E-4, 1.165028E-4, 1.171728E-4, 1.178466E-4, 1.185243E-4, 
    1.192059E-4, 1.198914E-4, 1.205809E-4, 1.212743E-4, 1.219717E-4, 
    1.226731E-4, 1.233786E-4, 1.240881E-4, 1.248017E-4, 1.255194E-4, 
    1.262412E-4, 1.269672E-4, 1.276973E-4, 1.284317E-4, 1.291702E-4, 
    1.299130E-4, 1.306601E-4, 1.314115E-4, 1.321672E-4, 1.329273E-4, 
    1.336917E-4, 1.344605E-4, 1.352337E-4, 1.360114E-4, 1.367936E-4, 
    1.375802E-4, 1.383714E-4, 1.391671E-4, 1.399674E-4, 1.407724E-4, 
    1.415819E-4, 1.423961E-4, 1.432149E-4, 1.440385E-4, 1.448668E-4, 
    1.456999E-4, 1.465378E-4, 1.473805E-4, 1.482280E-4, 1.490804E-4, 
    1.499377E-4, 1.508000E-4, 1.516672E-4, 1.525394E-4, 1.534166E-4, 
    1.542988E-4, 1.551861E-4, 1.560786E-4, 1.569761E-4, 1.578788E-4, 
    1.587867E-4, 1.596999E-4, 1.606183E-4, 1.615419E-4, 1.624709E-4, 
    1.634052E-4, 1.643449E-4, 1.652900E-4, 1.662405E-4, 1.674611E-4, 
    1.702246E-4, 1.730337E-4, 1.758892E-4, 1.787918E-4, 1.817423E-4, 
    1.847414E-4, 1.877901E-4, 1.908891E-4, 1.940392E-4, 1.972413E-4, 
    2.004963E-4, 2.038049E-4, 2.071682E-4, 2.105870E-4, 2.140621E-4, 
    2.175947E-4, 2.211855E-4, 2.248356E-4, 2.285459E-4, 2.323175E-4, 
    2.361512E-4, 2.400483E-4, 0.0, 0.0, 0.0};
    sbfill(sb, nbt, delta, be, ap, betan, ndmax );
    for (int i = 2429; i < 2679; ++i){
      REQUIRE( correct[i-2429] == Approx(sb[i]).epsilon(1e-6) );
    }
    for (int i = 0; i < 2429; ++i){
      REQUIRE(0 == Approx(sb[i]).epsilon(1e-6));}
    for (int i = 2679; i < ndmax; ++i){
      REQUIRE(0 == Approx(sb[i]).epsilon(1e-6));
    }
              
    for ( size_t i = 0; i < ap.size(); ++i ){ ap[i] = (i+1)*0.0001; }
    for ( size_t i = 0; i < sb.size(); ++i ){ sb[i] = 0.0; }

    be = 0.015;
    betan = {0.015, 0.018, 0.022};
    correct = {2.5442136E-4, 2.1217272E-4, 1.5127199E-4, 9.9999997E-5, 
      6.6106083E-5, 4.3700144E-5, 2.8888454E-5, 1.9097026E-5, 1.2624296E-5, 
      8.3454283E-6, 5.5168360E-6, 3.6469643E-6, 4.0461258E-6, 6.1097007E-6, 
      9.2257245E-6, 1.3930959E-5, 2.1035923E-5, 3.1764507E-5, 4.7964804E-5, 
      7.2427455E-5, 1.0936636E-4, 1.6514458E-4, 2.1798446E-4, 2.6092250E-4};
    sbfill(sb, nbt, delta, be, ap, betan, ndmax );
    for ( int i = 2541; i < 2565; ++i ){
      REQUIRE( correct[i-2541] == Approx(sb[i]).epsilon(1e-6) );
    }
  } // GIVEN
  GIVEN( "invalid inputs" ){
    WHEN( "not even close to enough entries in sb, ap" ){
      int ndmax = 100, nbt = 2545;
      double delta = 1.791435E-3, be = 0.15;
      std::vector<double> sb (ndmax,0), ap (ndmax,0), betan {0.15, 0.18, 0.22};
      for ( size_t i = 0; i < ap.size(); ++i ){ ap[i] = (i+1)*0.0001; }
      THEN( "an exception is thrown" ){
        REQUIRE_THROWS( sbfill(sb, nbt, delta, be, ap, betan, ndmax ) );
      } // THEN
    } // WHEN

    WHEN( "just barely not enough entries (1 too few)" ){
      int ndmax = 5088, nbt = 2545;
      double delta = 1.791435E-3, be = 0.15;
      std::vector<double> sb (ndmax,0), ap (ndmax,0), betan {0.15, 0.18, 0.22};
      for ( size_t i = 0; i < ap.size(); ++i ){ ap[i] = (i+1)*0.0001; }
      THEN( "an exception is thrown" ){
        REQUIRE_THROWS( sbfill(sb, nbt, delta, be, ap, betan, ndmax ) );
      } // THEN
    } // WHEN
  } // GIVEN
} // TEST CASE
